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We investigate the dynamics of a single phonon (oscillator) mode linearly coupled to an electronic 
few-level system in contact with external particle reservoirs (leads) . A stationary electronic current 
through the system generates non-trivial dynamical behaviour of the oscillator. Using Feynman- 
Vernon influence functional theory, we derive a Langevin equation for the oscillator trajectory that 
is non-perturbative in the system-leads coupling and from which we extract effective oscillator 
potentials and friction coefficients. For the two simplest cases of a single and two coupled electronic 
levels, we discuss various regimes of the oscillator dynamics. 

PACS numbers: 71.38.-k, 73. 21. La, 85.85.+j 



I. INTRODUCTION 

Nano-electromechanical systems (NEMS) are test- 
beds for the observation of fundamental quantum be- 
haviour of objects which are huge on the scale of indi- 
vidual atoms. For example, recent experiments^ have al- 
lowed a detailed study and control of single phonons by 
cooling a macroscopic resonator mode close to its ground 
state and coupling it to single electronic degrees of free- 
dom. 

One fascinating aspect of NEMS is their concep- 
tual simplicity that nevertheless can give rise to highly 
complex physics, and the links that can be estab- 
lished to other fields such as molecular electronics or 
optomechanics^i^. One of the challenges are the de- 
tails of the oscillator-electron coupling, i.e. in the lan- 
guage of measurement theory to understand, utilize^ and 
control^i^ the 'back-action' effects of the detector (e.g., 
superconducting single-electron transistors^!^) onto the 
oscillator. 

In many cases, even if no further approximations in 
simple theoretical models are made, the coupling to 
external electronic reservoirs that link the detector to 
the outer world is treated perturbatively, i.e., in the 
framework of (quantum) Master equations. This has 
turned out to be a highly successful approach, in par- 
ticular to describe such various systems as NEMS cou- 
pled to single electron transistors (SETS)^"— , Franck- 
Condon blockades!^ in transport through molecules 
with strong electron-phonon coupling, or quantum 
shuttlesi^"— . From the theory of electronic transport 
through nanostructuresii, however, it is known that 
such approximations usually are reliable only in the 
limit of high external voltage bias, where non-Markovian 
effectSiiS due to quantum coherences between the exter- 
nal reservoirs and the electronic system (SET, quantum 
dot etc.) can be neglected. It is therefore desirable to de- 
velop tools that allow a description of NEMS beyond the 
Master equation regime (weak electron-leads coupling) 
and at the same time are not merely perturbative in the 
coupling of the oscillator to the electronic environment^^. 

In the past, the coupling of electrons to a single bosonic 
mode has been solved exactly for the case where only one 



single electron is presen t-^ i.e. in an empty band ap- 
proximation. The inclusion of Fermi sea reservoirs at dif- 
ferent chemical potential transforms this into a difficult 
many-body problem out of equilibrium, and approxima- 
tions are necessary 
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Our approach in this paper is to combine exact solu- 
tions of the electronic system with a semi-classical ex- 
pansion, together with an adiabatic approximation for 
the oscillator dynamics within the Feynman- Vernon in- 
fluence functional (double path integral) theory^^. We 
revise this method, which has first been used for simple 
NEMS models by Mozyrsky and co-workers^^, and ex- 
tend it to allow for the description of a relatively large 
class of non-equilibrium electronic environments. The 
key idea is a systematic expansion around the classical 
path in order to obtain a Langevin equation for the oscil- 
lator. Already at the simplest level of this approximation 
(neglecting quadratic fluctuations around the diagonal 
path in the reduced density matrix of the oscillator), the 
coupling to the electronic non-equilibrium environment 
gives rise to non-trivial effects such as effective oscilla- 
tor potentials and non-linear friction coefficients lead- 
ing to both positive and negative damping21. Gaussian 
fluctuations around the classical path are built into the 
theoretical description, but they have to be evaluated by 
numerical solutions of the underlying Langevin equations 
which is not done in this paper. 

We compare two non-interacting electronic 'quantum 
dot' models with one and two levels between source and 
drain reservoirs: a single dot, and two dots in series. The 
oscillator couples linearly to the dot occupation (single 
dot) or to the occupation difference (double dot). One 
particular feature of the double dot case (where quan- 
tum superpositions of the electrons become important) 
is the occurence of limit cycles in phase space caused by 
a negative damping. 

The paper is organized as follows: after introducing 
the path integral formalism with a generic model in Sec. 
II, we present the single dot case in Sec. Ill and the 
double dot case in Sec. IV. Detailed derivations of the 
important formulae can be found in the appendices. 
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II. GENERIC MODEL 

A large class of NEMS can be described as a compo- 
sition of an electronic system "He, a mechanical system 
■Hose and a linear coupling between the two. Thus we set 
up a generic Hamiltonian "Hgcn by 



2m 



h -Hose - Fq, 



(1) 
(2) 



Here, Hose describes a single oscillator with p (momen- 
tum) and q (position) operators. The oscillator mode 
is confined in a potential Vosdq), F denotes an elec- 
tronic force operator, and m labels the oscillator mass. 
Whithin this paper the reduced Planck constant is set to 
one {h= 1). 

Our generic model does not include an additional os- 
cillator damping mechanism . In the usual Master equa- 
tion treatment of NEMS, Lindblad-form damping due to 
external degrees of freedoms is included phenomenologi- 
cally. In the path integral formalism used here, such de- 
grees of freedom can be easily included at least for linear 
or weak coupling to the oscillator. In order to elucidate 
the effect of the electronic environment that we treat in 
all orders in the coupling to external electronic reservoirs 
(contained in He), we choose not to include additional 
damping terms in our model here. 



introduced the angular oscillator frequency by and 
electron transition rates by Fl, we have to satisfy the 
condition cjq ^ Fl, Fr. 

In the next step, we derive a stochastic equation of 
motion for the classical trajectory, taking into account 
the propagation of the initial reduced density matrix, 
cf. appendix [X] The key step here is a cluster expansion 
to quadratic order in the off-diagonal path yt that de- 
scribes the Gaussian fluctuations around the classical os- 
cillator trajectory, where the fluctuations are determined 
by the properties of the non-equilibrum environment. To 
achieve a self-consistent equation of motion, we then re- 
insert the full time-dependence of the fixed classical tra- 
jectory in accordance with the adiabatic approximation 
and end up with the Langevin equation 



mxt + Ksc(a;t) - {F[x\{t))+xtA[x]{t) = 6- (5) 

Here the interaction picture of the electronic operator 
is given by F[x\{t) = exp[i(-Hc - Fx)t\F eyi-p[-i{'Hc - 
Fx)t\, and is a stochastic force with zero mean and 
the correlation function 



(66') -2Re(<5F[x](^)(5F[x](^'))• 



(6) 



The fluctuation of the electronic operator is defined by 
5F[x\{t) = F[x\{t) - {F[x\{t)). The friction A[x\ is given 

by 



A. Stochastic equation of motion 

We describe the oscillator dynamics by the reduced 
density matrix of the oscillator in position representation 
Posc{q,q' ,t) = {q\Posc{t)W) , for which we derive a semi- 
classical equation of motion for the oscillator position by 
using Feynman- Vernon influence functional theory sim- 
ilar to Mozyrsky and co-workers^. Assuming that the 
total density matrix x(t) factorises at the initial time 
into a system and a bath part x(to) = Posc(io) ® Pb, the 
propagation of the reduced oscillator density {q\posc{t)\q'} 
at time t is given by a double path integral3^^, cf. ap- 
pendix [Xj A transformation to center-of-mass and rela- 
tive coordinates 



qt + q't 



vt^qt- qt 



(3) 



A[x\{t)^2 dt' t' lm{6F[x\{t)6F[x\{t')). (7) 

Jto 

In the following we choose ti^ = —oo as initial time. For 
the specific cases of single and double dots, we checked 
that the upper integration boundary can be extended to 
infinity. 



III. ANDERSON-HOLSTEIN MODEL (AHM) 

The AHM combines a single bosonic mode with a sim- 
ple electronic transport system. We describe the bosonic 
part in first quantisation as 



He 



1 

2m 



1 



P 



2 -2 



(8) 



has the notion to detach the classical trajectory xt from 
the quantum mechanical deviations. Within a Born- 
Oppenheimer approximation, change of variables allows 
us to study a slow oscillator by an adiabatic approxima- 
tion of the classical trajectory 



Xt~XQ+ txo, 



(4) 



cf. appendix |^ In this approach, the typical timescale 
of the oscillator movement is slow compared to the elec- 
tronic transition rates. In the subsequent, when having 



with the bosonic position and momentum operators x 
and p, the phonon frequency loq and the phonon mass m. 
The oscillator length is defined by Iq = [muJoY^^- The 
electronic part is a single dot level confined between two 
leads: 



ka 



E 

ka. 



The dot level has energy and creation/annihilation 
operators d'^ /d. The operators cJ.^/cfeQ create/anihilate 



3 



C/eil X 10 



gr^/r 



Ml 


--^---x = 


occupation zero 
















-^Hooke p 




average occupation 






glo 




occupation unity 



FIG. 1: Sketch of the AHM model. Two leads (grey bars) with 
chemical potentials /iL, Mr = A*l — eVsias embed the dot level 
which couples to an oscillator. Here Vsias is the bias voltage 
and e the electron charge. The level energy Sd is shifted by 
the oscillator position x, the shifted level energy = Ed — Aa; 
is stationary, if the oscillator force Fnookc = mujgx and the 
electron force Fc = \{h[x]) are in balance. Here m and ojq 
describe the oscillator mass and angular frequency, and A (g) 
denotes the coupling constant, cf. eQ. pO|l . For sufficiently 
small tunneling rates Fl, Fr the stable stationary level en- 
ergies (dashed lines) are located at x = 0, ^ZoFl/F, grZo- The 
instable points of Sx (dotted lines) are located around the 
chemical potentials. 



electrons with momentum k and energy Ska in a free elec- 
tron gas in the left {a = L) or right (a = R) lead. Elec- 
tronic transitions are possible with amplitude Vka be- 
tween the dot and a state in the lead. Between the two 
subsystems there is a simple linear coupling with coupling 
constant A, such that the total Hamiltonian reads 
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FIG. 2: left) Density plot of the effective oscillator poten- 
tial (7eff as a function of oscillator position (a;) and Vsias 
in units of and Iq with the parameters Fl^/i^o = 0.7, 
£d/i^o = 2.9, Puq = 10 and g = 2.4. For increasing bias 
voltage, the effective potential shows two, three, two and one 
minima. This different regions are separated by white lines 
at VBias/t^o = 2.04/2.58/2.64. Right) Uc« and the two con- 
tributions -Fkooko, to the force F^a at bias voltages (sym- 
metric choice) Vsias/i^o = 0.0/2.3/5.0. UcS exhibits extrema 
where the two force contributions (grey and dashed line) are 
in balance. The width of the center plateau in the electronic 
force contribution Fc a (n) grows with increasing bias volt- 
age Vsias, cf. Fig. Il]) for the intermediate occupation (n). For 
sufficiently high bias voltage only one minimum remains. 



"Hahm = T-Lc + "Hose - Xci^di 



(9) with the effective force Fcs and the friction A[x\: 



For convenience we introduce the dimensionless coupling 
constant: 



A 



mujQlo 



(10) 



Here, we regard spinless electrons. Note that the general- 
ization to a model including both spin directions requires 
an onsite interaction term like Ufi^fi^. In the following, 
we only consider non-interacting electrons. A physical 
realization of the model discussed here would correspond 
to either spin polarized electrons or a coupling to the 
oscillator that effects only electrons of one certain spin 
direction. 



A. Langevin equation 

When applying our generic form ^ to the AHM, we 
obtain the Langevin equation 



Fcs{xt) = -mujf^xt + X{n{t)) = — — 
A[x]{t) = 2A2 f dt' t'lra{Sn{t)Sn{t')). 

J to 



(12) 



The effective force has two contributions; a term propor- 
tional to the elongation (Hooke's law) and the electron 
force which is proportional to the dot occupation; which 
only contributes if the dot is occupied. The friction term 
a;tA[x] results from stochastic electron jumps between the 
leads and the dot. 

For finite bias voltage, figure [T] shows the positions of the 
dot energy level and the points of instable balance, which 
result from the balance of both forces. 

The effective force and therewith the oscillator poten- 
tial are determined by the dot occupation n(t), whereas 
the friction and the stochastic force correlation {£,t£,t') = 
X'^2Re{Sh{t)Sh{t')) depend on the imaginary /real part 
of the dot correlation function. The dot correlation func- 
tion in terms of the lesser and greater Green's function 
(which are derived in appendix |B]) reads 



mxt - Fes{xt) + XtA[x]{t) = 



(11) 



{6n{t)Sn{0)) = G<{~t)G>(t). 



(13) 
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B. Effective potential 

The occupation of the dot is calculated in an adia- 
batic approach with the help of the lesser Green's func- 
tion G^{uj), cf. appendix iBl where we assume constant 
tunneling rates Fq,. For finite temperatures, the dot oc- 
cupation reads 

= -i— ( duj G<{uj) 



> — Im^'- + — 



1 /3F . /3(£^-Ma) 
h i 

2 47r 27r 



(14) 



whereby /3 stands for the inverse temperature, Ex is a 
short notation for Ed — Xx, fj,a denotes the chemical pote- 
tials and F = Fl + Fr . 5" designates the Digamma func- 
tion. By integration we obtain the effective potential 



C/eff(x) = - / dx' F,s{x') = 
Jo 

/3(£d - Xx - 



(15) 



lnF|^^ + i 

/3(£d - Ma) 

2tt 



271 



1 X 

2k 



X 



Wo, 

(16) 



with F(-) denoting the Gamma function. In the absence 
of friction the dynamics of the oscillator is determined 
by the effective potential Uas{x). The high temperature 
case is of minor interest, because the temperature washes 
out the structures of Ucs{x). 

In Fig. [2] we present features of the effective potential 
Ucs(x) at zero temperature. Ues(x) has minima when 
the effective force Fee{x) is zero. This is the case when 
the oscillator force and the electron force are in balance. 
We plot the two contributions to the effective force: the 
force ^Hookc is proportional to the displacement (x) , the 
contribution results from the dot occupation (scaled 
with g) and has two steps at 



1^ 



£d 

Wo 



Wo 



a e {L, R}. 



(17) 



In the upper part of Fig. [5] we easily recognize that for 
large bias voltage there will be only one minimum at 
x/li) = ^Fl/F. By changing the bias voltage or the cou- 
pling strength, we can reach situations where two minima 
at x/Iq = and at x/l^ = g are added, like in the middle 
part of Fig. [2] and also situations, where only the two 
minima at the side remain and the one in the middle 
vanishes, like in the bottom part of the [5] (this agrees 
with Mozyrsky et ali^). Increasing bias voltage shifts 
the steps in F^ apart, whereas increasing the coupling 
constant minimises the distance. The positions of the 
minima/steps are exact in the zero rate limit (Fl — > 0, 
Fr — )■ 0) where the averaged occupation is step-like. For 
finite rates, the steps smoothes out. 



The form of the oscillator potential already hints to- 
wards the phase space spanned by position {x) and mo- 
mentum m{x). In the absence of friction, the minima of 
the potential correspond to the fixed points of the oscil- 
lator motion. 



C. Phase space portrait 

The classical trajectory {xt) is obtained by neglecting 
fluctuations due to the stochastic force in pT|): there- 
fore the equation of motion reads 

m{it) + {xt)A[{x)]{t) - F,s{{xt)) = 0. (18) 

In Fig. [3] we show the effective potential and the fric- 
tion, as well as the solutions of the classical equation of 
motion in phase space without and with friction. The 
rows correspond to three different voltages leading to the 
three typical cases with one, three and two minima of 
the potential (the latter case is investigated in^i). The 
initial conditions are chosen for each phase diagram sep- 
arately in order to make the characteristic shapes visible. 
The phase diagrams without friction follow directly from 
the shape of the effective potential. The trajectories in- 
cluding the friction follow the ones without friction for 
a while until a position with a peak in the friction is 
reached that produces a kink-like damping feature. The 
friction causing the kinks is maximal, when the shifted 
energy level Ex is in resonance with the chemical poten- 
tials ml,r — ±^Bias/2 (Fig. [T]), because at these points 
the average occupation switches and the electronic fluc- 
tuations and therewith the friction itself is very large. 

Between its peaks the damping does not completely 
vanish at finite bias, so all phase space trajectories end up 
in spirals and reach stable fixed points after infinite time 
(we have terminated all phase trajectories after the same 
time). The position of the peaks (instable fixed points) 
separates the stable fixed points. In the two minima case 
(^Bias/i^o = 0.0 in Fig.[3l) the left fixed point corresponds 
to the zero occupied electronic level in Fig.[T]and the right 
fixed point to the fully occupied level. For increasing bias 
voltage, the friction splits and we obtain a third fixed 
point in between, corresponding to the average occupied 
state. For sufhciently high bias voltage (large transport 
window) only the average occupied state survives. 

Apart from the stable fixed points we also observe 
saddle points that repell the trajectories near the kinks 
where the damping is large. They result from the kinks 
in the dot occupation (Fig. [21 right) and correspond to 
the dotted lines in Fig. [T] and the bright V-structure in 
Fig. [51 left. At large bias the saddle points move out of 
the range of allowed positions, the same happens to the 
peaks in the friction. At infinite bias the effective poten- 
tial becomes a simple parabola and the friction vanishes 
completely (this can be checked analytically). 
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FIG. 3: Effective oscillator potentials (A), frictions A{x) (B) and the phase space portraits without (C) and with (D) friction 
in units of ojq and lo with the parameters Fl.r/ijJo = 0.7, £d/wo = 3.0 and g — 2.45 at zero temperature. From the bottom to 
the top the values of the bias voltage (symmetric choice) read yBias/<^o ~ 0.0/2.2/6.0 in correspondence to the three regions in 
Fig. [2] (Right). The peaks of the friction are located at j: = [sd T VBias/2]/A where the shifted level energies Ex are in resonance 
with the chemical potentials fJ.L,R, cf. Fig. [T] 



IV. DOUBLE QUANTUM DOT SYSTEM (DQD) 

The DQD that we treat in this section consits of two 
single dot levels coupled by a tunnel barrier. Again we 
assume a coupling to a single bosonic mode. The total 
Hamiltonian is composed of the oscillator part i?osc, cf. 
eq. ([8]), the electronic part and an interaction part 
which describes the coupling between the oscillator and 
the two dots. In contrast to the AHM, here the oscillator 
couples to the difference of the occupation numbers with 
the coupling strength A. The total Hamiltonian therefore 
reads 

Hbqd = i?c + Hose ~ Xx{dldi^ ~ 44), (19) 
containing the electronic part 

ka. a 

+ E [^k<^^^lJo' + yLdUua] ■ (20) 

ka 



In this model, d\/da denotes the creation/annihilation 
operator of the ath [a e {L,R}) dot. designates the 
corresponding level energy and Tc describes the tunnel 
coupling matrix element between the two dots. Note that 
we include no Coulomb interaction terms here. 

The stationary average of the population difference 
(ctz) corresponds to the electronic force operator F of 
the generic model, eq. ([T]). The occupation of the ath 
dot 



{na) ^ I dco G<Jio) (21) 



can be calculated by using Keldysh's equation (jC5l) . see 
section IC II Therewith the occupation difference of the 
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FIG. 4: Examination of the DQD system for various tunnel couplings Tc, increasing from left to right. Row A shows the results 
of a fixed point analysis. A corresponds to the determinant and r to the trace of the Jacobian matrix. Rows B and C display 
phace space portraits without and with friction. Explicit parameters are iTcl^ = 0.4; 0.43; 0.49; 1.0 Wq. We used equal rates 
Fl = Fr = l.ScjQ. For the chemical potentials we assumed /xl = Ihojo and /ir = —5hjJo. The dimensionless coupling constant 
is chosen as g — 2.5, and the internal bias voltage as Vint = Bhuo, wheras i/l = — = eVintfi. 



left /right dot follows from 



leads to 



F 

47r 



duj 



w4 + 2Aa;2 + 



,(22) 



with the abbreviations 14,3 — Vh.R T Aa; and 



A = - [|r,p + t>g - (r/4)2] , s = [|T,p + + (r/4)2] , 

(23) 



r 

47r 



E 



2J^L2:2(/ia,A B) 



+ sgn(z/QZ^L)|2:3(Ma,^,5) + 

(?^ + (r/4)2-|Tcp)Xi(Mo,AB) 



(24) 



They are expressed in terms of the auxiliary functions Tj 
defined in appendix iDl 



A. Langevin equation 

When applying the generic model, cf. eq. ([5|), to the 
DQD Hamiltonian -ffogo, we obtain the Langevin equa- 
tion 



whereas we assumed — —Vh- Calculating the integrals 



mxt - Fesixt) + XtA[x]{t) = 



(25) 
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with the effective force Fcs and the friction In 
contrast to the AHM, the effective force is affected by 
the population difference and not only by the occupation 
number. Explicitly, 



(26) 



A[x]{t)^2X^ dt' t'lia{5d,{t)5d,{t')). (27) 

J to 



In appendix I C 2 1 we derive the explicit expression for the 
friction. In the case of infinite bias (IB) and ?r = — J'l 
we obtain 



^iBN(t) --8 iT^r — 



(28) 



In contrast to the AHM the friction does not disappear 
for infinite bias. The second difference that we can recog- 
nise by regarding the prefactor i7l is that we obtain re- 
gions where the friction is negative. The latter also holds 
for the finite bias case. 

The real part of the correlation function according to 
the population difference az = dl^di^ — dpdn determines 
the correlation function of the stochastic force 

(66') = X^2Re{6az[x]{t)Saz[x]{t')), (29) 

with fluctuation (Jti^ [x] (t) = (T2[a;](t) — {az[x]{t)). 

B. Fixed point analysis 

The effective potential determines the behavior of the 
oscillator trajectories in the phase space ((p)-(a;)-plane), 
as seen in section IIIIBI for the Anderson Holstein model. 
In the following we examine the differential equation of 
the system by studying its fixed points^i^. Some fur- 
ther theoretical details of this kind of investigation are 
explained in appendix IC 31 

For the double dot we obtain the dynamical system 



(Pt) 



-iPt) 
m 



lo 



lo 



.9(^40) 



iiA[{xm 



•(30) 



Fixed points occur under the condition (pt) = (xt) = 0, 
i.e. (pt) = and following from that, the fixed points 
position coordinates are equal to the roots of the effective 
force 







+ g{az{t))^0. 



The Jacobian matrix is obtained from 

/ 1 



J* = 



"0 

lo 



evaluated at the fixed point (x*), whereby {p* 
Determinant A and trace r become 



A = — 



UJQ 

lo 



1 

1- ~ .9 



d 



d{xt) 



-A\ 



(31) 

The trace decides about the stability of a fixed point 
and is equal to the negative friction here. For the case 
without friction the trace is equal to zero, therefore only 
centers occur in the phase plane. In the case with friction 
the trace can be either positive or negative leading to 
both, stable and unstable fixed points. For comparison 
see Figure SI where row A depicts the results of a fixed 
point analysis. 

The effective force -Fcff is plotted in the upper part 
of each plot in row A together with trace and determi- 
nant. Therewith the characteristics of the fixed points 
are determined and it is possible to predict the shape of 
the phase space portrait. In the lower parts of the plots 
in row A these predictions are illustrated, fixed points 
are marked by black lines. For small tunnel coupling 
(diagram Al) we obtain seven fixed points. These can 
be characterised as three stable and one unstable spirals, 
each seperated by one of the three saddle points. Sta- 
ble spirals correspond to the different rest positions for 
the oscillator and the saddle points to the points of in- 
stable balance. Increasing the tunnel coupling leads to 
five fixed points. In graph A2 we obtain a stable spiral 
in the middle enclosed by two saddle points and followed 
by a stable spiral on each side. In the diagrams A3 and 
A4 the mean point changes to an unstable spiral. With 
further increased tunnel coupling (A3) there remain only 
three fixed points and for even higher values of Tc (A4) 
the number of zeros in the effective force reduces to one. 
For the oscillator the latter means that it is shifted to a 
new rest position independent from its inital position. 



C. Phase space portraits 

In Figure 13] the rows B and C depict phase space por- 
traits for the double dot system without and with fric- 
tion. The four columns correspond to four different val- 
ues of the tunnel coupling Tc- As initial condition, the 
momentum was set to zero and the positions were chosen 
in order to show the various shapes of the trajectories. 
In the case without friction we recognize periodic cycles 
which are stable and run around one or more fixed points. 
These centers were also expected from the analysis in sec- 
tion |IVB] The fixed points correspond to certain states 
of the double dot system. If the left dot is occupied 
the rest position of the oscillator is shifted to the right 
and correspondingly to the left for an occupied right dot. 
These points turn to stable spirals when the friction is 
turned on. The states when both dots are occupied or 
empty correspond to the fixed point in the middle. 

The lowermost row shows what happens when we in- 
clude the friction in our calculations. In contrast to the 
single dot system, here the friction has positive as well 
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as negative values depending on the position of the os- 
cillator. This means that the oscillator is either deceler- 
ated or accelerated. Both can be interpreted as inelastic 
jumps of the electrons, where energy is transfered be- 
tween the electrons and the oscillator in both directions 
like it has been observed in^i. There, the authors consider 
a resonator coupled to a superconducting single electron 
transistor (SSET). As a result of the interplay of positive 
and negative damping in a certain parameter range they 
observed limit cycles and bistability in the phase plane, 
in our work we obtain a likewise behaviour for the os- 
cillator. In contrast to our semiclassical approximation 
they investigate the Wigner function of the system with 
numerical master equations. They compare these results 
to a mean field evaluation of the expectation value of 
the oscillator position'^^. For weak coupling the mean 
field approach gives quantitatively correct results, and 
for higher coupling it still describes the dynamics qual- 
itatively correct. These results suggest that our use of 
average oscillator positions and momenta is qualitatively 
correct for a description of the oscillator instabilities. 

Consider again row C in Figure lU where the results 
for the dynamics of the oscillator with friction are plot- 
ted. The outer left and right stable spiral do not change 
by increasing the tunnel coupling Tc- By contrast the 
oscillator's behaviour between these stable rest positions 
changes a lot. In graph CI we observe a stable and an 
unstable spiral, like we expect from the fixed point anal- 
ysis. In the neighbourhood of the unstable fixed point 
the friction is negative, so the oscillator trajectory is re- 
pelled and ends up in the right stable spiral. In diagram 
C2 we recognize that the latter path becomes stable. The 
limit cycle appears when the unstable spiral has disap- 
peared and exists as long as the tunneling coupling is 
in the range of 0.42 < iTdV^o < 0.5. There we ob- 
serve a bistability: as the initial position gets closer to 
the fixed point the limit cycle turns into a stable spi- 
ral. By further increasing Tc, the middle spiral becomes 
unstable and a second limit cycle appears (C3). This 
limit cycle with a smaller radius exists in the range of 
0.49 < iTclVw^ < 0.54. For |Tc|/wg ~ 0.49 the system 
undergoes a Hopf bifurcation^^, which happens when a 
pair of complex eigenvalues from the dynamical system, 
which determine the evolution in the phase plane, see 
section fC 3[ cross the imaginary axis from the left to the 
right half-plane. In other words, the trace t changes 
its sign and at the bifurcation point the eigenvalues are 
purely imaginary Ai,2 — =F*2\/A, see equation (jC16|) . In 
the last graph of row C both limit cycles have disappeared 
and the oscillators path ends up in the left or right stable 
point. 

In our calculations we choose F = Swo, standing in 
some contrast to our adiabatic approach, which implies 
a slow oscillator (F ^ wq). If F ^ wq the interaction 
between the current and the oscillator is strongest be- 
cause both act on the same timescale. Interesting effects 
still appear with a slightly enlarged F like in our plots, 
but for F 3> Wo there remains only one stable spiral. This 



means, that the oscillator rest position is shifted from its 
ground position caused by the stochastic processes initi- 
ated by the current. We presume that our approach is 
useful also for a comparative fast oscillator and we will 
accomplish further investigations with a non-adiabatic 
approximation to reconsider our results. 



V. SUMMARY 

We have derived a stochastic equation of motion that 
describes the dynamics of a single oscillator coupled to 
an electronic environment out of equilibrium. We studied 
two cases, namely the single dot level and double dot 
(two-level) electronic system. For both cases we have 
explained the features of effective potential and friction 
for the ensemble averaged oscillator motion. The effects 
we recognize fit well together with former works. In the 
DQD model limit cycles and bistabilities appear. 

Until now the master equation has been used for most 
investigations of the oscillator behaviour in NEMS. We 
have used a method that gives us access to regions where 
the master equation has problems: we naturally include 
finite bias, and arbitrary electron coupling to external 
reservoirs. 

We had to stay in a regime with a relatively fast oscil- 
lator in order not to miss the interesting physical effects. 
The validity of our method in this regime could still be 
improved with a non-adiabatic calculation. 
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Appendix A: Path integral representation of the 
reduced density matrix 

The reduced density matrix p{t) = trBx(i), obtained 
from the total density matrix x{t) by tracing out the 
bath degrees of freedom, describes the mechanic sub- 
system "Hose, cf. ([T]). By using a factorising initial 
condition x(^o) = p{to) Pb the elements of the re- 
duced density matrix can be expressed in a path integral 
representation22i^ 

(glPosc(Ok') = J dqo dq'o {qo\Poscito)\q'o) 

r<iH) r<i'{t) 
X Vq{T) V*q'iT) e'^^^~^^'^T[q,q']iT). 

Jq{to) Jq'ito) 

(Al) 

Hereby Sq = j^^ dt' [-^mq^, — Ksc(<7t')] denotes the clas- 
sical action in the path integral for eq. (jAip and the 
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Feynman- Vernon influence functional is defined by 

F[q, <z'] (to , t) - ti-B { [q] {t, to)U[q] (t, to)pB } • (A2) 
The time-evolution operators are defined by 

(7[g](t,to) = T^exp[-z f dt' HreskKt')] (A3) 

Jto 

corresponding to the reservoirs 'Hrcs[q]{t) = "He — Fqt- In 
the next step we establish an effective interaction picture 
by 

U[q]it,to) ^ e^'^"^'-'"^ U[q]{t,to) 



-T^cxp[-z / dt' V[q]{t')], 
Jto 



_ i-Ho{t-ta) 



V[q]{t)e 



-i-Hoit-to) 



(A4) 



Hereby the reservoir is decomposed in an unperturbed 
part Ho — Ti-c — Fxq and the perturbation y[g](t) = 
—F{tio + \yt) and V[q']{t) — ~F{txo — \yt), where we 
wrote w xq + tx^ using the adiabatic approximation 

©■ 

By using eq. (jA4p . the influence functional can be ex- 
panded to second order in perturbation theory. Denoting 
( • > = ti-BlPB • }, Vi V[q']{ti) and Vi := V[q]{ti), the 
Feynman- Vernon influence functional reads 

T[q,q']{h,t) ^ {U[q']{tM)UmM)) 



^l + i I dh {Vi - Vi) 



+ / dh dh {{Vi - Vi)V2 - Vi{Vi - Fi)). (A5) 



When plugging in the definition of the perturbation (in- 
teraction picture) and taking into account that the force 
operator F is hermitian, we obtain 



F[x + y/2,x~y/2]{to,t) = l + i dh yt, 

Jto 

{F{h)} - 2io / ' dt2 t2lm{F{ti)F{t2)) 



t rti 

dh dhyuRe{F{h)F{t2))yt,. (A6) 

to Jto 

Performing a cluster expansionSS., the influence func- 
tional can be expressed in terms of the influence phase 



<^[x,y]{to,t) 



\nT[x + y/2,x^y/2]{to,t) ^ ~i / dh yt, 

Jto 

{F{h)) -2x0 dt2 t2lm{dF{ti)6F{t2)) 



f dh r dt2 yt, Re{SF{h)SF{t2))yt,, (A7) 

'to Jto 



whereby SF{t) — F{t) — {F{t)) denotes the fluctuation 
around {F{t)). This equation can be easily verified by ex- 
panding the exponential exp[— $[a::, y]{to, t)] to second or- 
der in perturbation theory and comparing with eq. (jA6l) . 

Together with the classical action, we define an effec- 
tive action functional by 

A[x,y]{to,t) := Sx+y/2 - Sx^y/2 + i^x,y]{to,t) 
pt 

= / dh [mxii/i - Ksc(a;i + ^) + Ksc(a;i - -^)] 
Jto ^ ^ 



+ i^\x,y\{to,t). 



(A8) 



whereby the reduced density matrix expressed in terms 
of the new variables, eq. ([3]), reads 

{x + ||p(i)|x '^^^ y lp(^o)|2^o - y) 



r^{t) py{t) 

/ Vx{t) / V*y{T) e 
Jx(to) Jy(to) 



JA[x,y](to,t) 



(A9) 



In a semiclassical approach we assume small deviations 
of the off-diagonal trajectories y from the diagonal ones. 
Thus the potential difference leads in second order in y 
to 

Voscix, + |) - l/osc(a:i - |) = VU^) + 0{y% 

(AlO) 

This approximation is exact for quadratic potentials, as 
we treat in this paper. Furthermore, with the boundary 
conditions y{to) = y{t) = and integration by parts the 
effective action functional is quadratic in y. 



A[x,y]{to,t) - / dti yi 



mi, + V^,,{x,)-{F{h)) 



+ 2io / dt2 t2lm{SF{h)SF{t2)) 

Jto 

+ 1 I dh r dt2 yiRe{6F{h)SF{t2)}y2] 

Jto Jto 



dh yi^iN 



to 



2./t, 



dh / dt2 yiLi,2[x\y2- 



(All) 



Completing the square of the integral kernel^^ 
^ D*yeyi'p{iA[x,y\}, the resulting path integral de- 
scribes a stochastic process with Langevin equation 

Kt[x] = mit + V:,,,{xt)-{F{t)) 

+ 2x0 I dt' t' lm{SF{t)SF{t')) = (A12) 

Jto 

with <^t a Gaussian stochastic force. To obtain a selfcon- 
sistent equation of motion, we finally replace Xo by Xt- 
This substitution concerns also the unperturbed Hamil- 
tonian which leads to Ho He — Fxt- The Langevin 



10 



equation then reads 

+ 2xt dt' t' lm{dF{t)SF{t')) ^ £,t. (A13) 

J to 

The term quadratic in y in eq. (jAlip determines the cor- 
relation function of the stochastic force, 

i^S') ^2Re{dF{t)dF{t')). (A14) 



For finite temperature one has to regard Fermi func- 
tions instead of the Heaviside theta function. By virtue 
of the residue theorem one finds 

2 TT^ r I2 47r ^ 27r j' 

(B6) 

with the digamma function 5'. 



Appendix B: Green's functions of the single dot 



The single dot lesser Green's function without coupling 
in energy space is^4 



(Bl) 



Here we have used the adiabatic approximation for the 
center of mass coordinate, thus the coupling to the os- 
cillator simply shifts the level £d by X{x), so we have to 
replace Ed by the shifted energy Ex = Sd — X{x). The 
greater Green's function is correspondingly 



.rL[i-/LH] + rR[i-/R(w)] 



(a;-£)2+r2/4 
F 



(B2) 



Hereby /a(w) = /(w - = l/[exp{/3(a; - Ha)} + 1] de- 
note Fermi functions with lead index a, inverse temper- 
ature /? and chemical potential fXa- The time-dependent 
Green's functions are obtained by Fourier transformation 
as G^{t) = f ^e~'''^*G^{uj). In the zero-temperature 
limit we have to replace fa{^) = ©(/^a — Then for 
t 7^ an integration leads to 



G^it) = ' 



27r ^ r 

a 

e+r/2|*l El { [tn^ sgn(0 + F/2] \t\} sgn(t) 
e-r/2|*l Ei{[in^ sgn(t) - r/2]|t|} sgn(i) 



±27rie"^/^l*le(±rJc 



where the first exponential integral is defined by 

^ — xt 



Ei(a;) 



dt 



'1 t 

For the case ^ = the lesser Greens function reads 



(B3) 



(B4) 



1 f 

,G<{t^O)^— di^Y. 



(tj-e)2-HF2/4 



(B5) 



Appendix C: Characteristics and calculations for the 
double dot 

1. Green's function of the DQD 

The Green's functions in the frequency domain accord- 
ing to the Hamiltonian in equation (|19p is derived via the 
equation of motion method. There the Green's function 
G is defined as resolvent of the Hamiltonian TJq via 



(wl - Hq)G{lu) 



a 



(Cl) 



Denoting the bath states with \(f)x) = |A) the previous 
equation yields 



(C2) 



Taking the matrix elements and inserting the Hamilto- 
nian leads to a set of equations, from whom the Green's 
functions are derived. 

Here, the electron-phonon coupling is described adiabat- 
ically, so the interaction part Hsb came in by shifting 
the dot level energies ^l.r J'l.r with i/l,r = j^l.rTAx. 
Finally the dot Green's function is given through 



Gd(^) 
with the elements 
Gll(c^) 



Gll(^) Glr(c^) 
Grl(w) GRR(a;) 

w - j7r - ^r{uj) 



[uj -ui^- SlCw)] [w - t/R - SR(a;)] - \Tc\ 

- i^L - - i^R - i;R(a;)] - \Tc\^ 

Glr(w) = — " G^r{uj) 



GRR(a;) = 



Grl(w) 



£c 

CJ - /?R - Sr(w) 



GLL(t^) 



(C3) 



Hereby we introduced the self energy Sq, (a € R, L) cor- 
responding to the left or the right dot with 



E„(a;) 



k 



(uj). 



(C4) 



9ka,ka IS the undisturbcd Green's function for the leads. 
We derive the associated retarded and the advanced 
Green's function by using the continuation rules 

G(w±iO+) ^ G^^^(a;) 
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The derivation of the lesser/greater Green's function is 
taking usage of the Keldysh equation: 

Gf, = E Gf;,H sfH G^^^H (C5) 

7 

whereas we assume both dots initially unoccupied. The 
lesser/greater self energy follows from 

k 

S>H = -zT^iuj) [1-AH]. (C6) 

2. Calculation of the friction for the DQD 

The friction is determined by the imaginary part of the 
correlation function of the double dot. From equation 
P7p. we start with expressing the correlation function 
through Green's functions and use their Fourier trans- 
forms: 

A[x]{t) 

= 2A2 [2'5a,^ - 1] / dt' t' 1 * 

[G<Ji' ~ t)G>,{t t') G<^{t - t')Gl^{t' t) 
= ^Im^ [25c,,/3 - 1] j duJi j duj2G'^^{i^i)* 

aui2 J 27r 

= ^ E [2'5-.^ - 1] / ^^G'<J^)|-G>,(^). (C7) 



In the last step a term was identified as the derivative of 
Dirac's delta, this result agrees with the solution in the 
work by Mozyrsky et al.— . 

The lesser/greater Green's function Keldysh equation fol- 
low from the Keldysh equation (|C5P and we consider the 
tunneling rates to be frequency independent, therefore 
the self energies are (T = 0): 

S>(t^) = <d{uj~ ^1^). (C8) 

Because of the lengthy calculation, we just outline the 
calculation for the LL-term Tll- The other terms 
(Trl, Tlr and 7rr) can be derived in a similar way. Here, 
the product of retarded and advanced Green's functions 
equates the squared modulus of Ga.p, so we can write 



G«lM G^l(^) = IGfj^)^ = \Gti.{u:)\\ 
G^u;) G^l(^) ^ \G^nMf = IG^J^)!^. (C9) 



For simplicity in the following we omit the superscript of 
the retarded Greens function and abbreviate the moduli 

by 



We obtain for the first term of the friction 



Tll = 



dio 



G^l(^)£g>l(c.) 



did 



^IGllHP 



|GLL(^)Pe(AiL - c^) + |GLR(c<.)pe(MR - u:) 
Q{u - ^l) + \Gi^l{uj)\^S{uj - //l) + 



^IGlrHP 



- /^R)|GLR(w)p5(a; - ^r) 



r 



(CIO) 



Some parts of the integration terms can directly be eval- 
uated with the help of the delta and Heavyside functions. 
By choosing the condition ^l > Mr, we get rid of a case 
distinction, which would be necessary for two product 
terms which included a delta-function. After performing 



also an integration by parts, we arrive at 

r2 



Tl 



ll 



|GLL(AiL)|' + |GLR(m)l' 



2 |Gll(a*l)|' |Glr(a*l)|' +4 |GLL(/iR)|' IGlrIa^r)!' 



4 J dcj|GLL(w) 



I^IGlrM^' 



(Gil) 
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In the following the Green's functions, derived in sec- 
tion [Ul] were inserted, whereas we assume equal tunnel- 
ing rates for the left and the right side, Tl = = ^F. 
Then a number of integrations by parts is performed to 
dispose of the derivation in the integral term. So we ob- 
tain a closed expression, whereas N{uj) abbreviates the 
denominator of the Green's function 



In an analogue way the other terms in (jC7|) can be de- 
rived, and finally the solution for the friction is 



r 



2 IT, 



16 



^(mr) 



^(mr)2 



PL 



MR 



7V(w)2 

(C12) 



A[x]it) 



TT 8 



-2 IT, 



-t- 



IT, 



+ 2 IT, 



4 IT, 



(a*r - '^r)(mr - J^l) - fe 



^(mr)^ 



8 |Tc|^ (i^L - i^r) I duj 



N{u;y 



(C13) 



r 



With the Assumption vr, = —I'h we get: 
N{uj) 

with A = -f |T,|' - rVi6), 
B= ?2 _^|T,|' + rVi6, 



(014) 



and the integral in equation (jC13p is given trough 
Ti{A,B) in section iDl 



3. Fixed point analysis for the double dot 

The fixed points of a nonlinear two dimensional system 
can be investigated with standard methods for linear dy- 
namical systent^. 

The general solution for a two dimensional linear system 
X = A X is 



x(t) = cie^i*vi C2e^='v2 



(015) 



and so determined by the eigenvalues \i^2 of the matrix 
A. The constants ci^2 depend on the initial conditions 
and Vi_2 are eigenvectors. 
The eigenvalues can be obtained from 



Ai,2 = i(T± Vt2-4A), 



(016) 



thereby r corresponds to the trace and A to the deter- 
minant of A. These two qualities determine the evolution 
of the trajectories in the phase plane. 
For a fixed point x* the condition x = must be ful- 
filled. Various classes of fixed points exist, whereas the 
determinant A decides which kind of point appears. In 
case of saddle points the determinant is negative and it 
is positive for spirals or nodes. The difference between a 
spiral and a node is that for the second one the eigenval- 
ues have no imaginary part. 

The trace r defines the stability of nodes and spirals, 
this is caused by the fact that t determines the sign of 
the eigenvalue's real part, for instance with negative real 
part decaying oscillations occur and the fixed point is sta- 
ble, see equation (j015|) . There also exist some borderline 
cases, whereas the centers are the most significant ones, 
they occur when the trace is equal to zero. 

This analysis can be assigned to a two dimensional non- 
linear system x = /(x). By assuming a small disturbance 
u = X — X* from a fixed point, we can invest if this distur- 
bance grows or decays by performing a Taylor expansion 



Ml 



U2 



fl[Xi,X2) +ui—\ -f M2— 4-h.t. 

OXi ^1 0X2 2 

j2{Xi,X2) +Ui—\ -HM2— +h.t. 

OXi '^1 0X2 2 



dfi 



(017) 
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The first term is zero and higher terms (h.t.) can be ne- 
glected becaiise the disturbance is small. So we get a 
linearised system u = J*u, containing the Jacobi matrix 
J* evaluated at the fixed point coordinates. The above 



explained analysis can be performed for this system. This 
is valid as long no borderline cases occur, then the higher 
terms may be more important. 



Appendix D: Auxiliary integrals (DQD) 



Ii{li,A,B) 



f 

J — oc 



duj- 



1 



1 



1 



1 



arctan ( 
arctan (- 



\/A + y/W^^ WA + y/W^W 



A - - B2 V ^ + \/A2 - B'^ 



du!- 



^ + 2Aoj'^ + 52 
1 



2VA2 - £2 



- i In (^2 _^ ^ _ ^^2 _ B-2) 
-^\n{n^ + A+ VA^ - B2) 



i{fx,A,B)= r 

J — cc 



duj- 



;4 + 2ylw2 + 52 



2VA2 - 52 



A - \/ A'^ - B'^ arctan | 



+ y''^4H-a/^^^^ arctan ( 



VA-\/l23Bf 

N 

\/a + 7^^b^' 
1 



) 



Ti{A,B) 



du)- 



j'^ + 2Auj^ + B^Y 



SB'^iB'^ - A2) 



2a;(52 _ 2^2 _ ^^2^ 

52 + 2Aa;2 + 
(A2 - 352 + y4^A2 - B^) arctan 



VA2 _ _B2 712^:^ 

(^2 - 352 _ A^/A^ - S2) arctan 



+ 



V^2 _ 52 + VA2 - S2 



(Dl) 



14 



* |robe rt@itp.physik.tu-berlin.del 

^ A. D. O'Connell, M. Hofheinz, M. Ansmann, R. C. Bial- 

czak, M. Lenander, E. Lucero, M. Neeley, D. Sank, H. 

Wang, M. Weides, J. Wenner, J. M. Martinis and A. N. 

Cleland, Nature 464, 697 (2010). 
^ S. Groblacher, J. B. Hertzberg, M. R. Vanner, G. D. Cole, 

S. Gigan, K. C. Schwab, and M. Aspelmeyer, Nat. Phys. 

5, 485 (2009). 

^ F. Marquardt and S. M. Girvin, Physics 2, 40 (2009). 

* A. Naik, O. Bun, M. D. LaHaye, A. D. Armour, A. A. 
Glerk, M. P. Blencowe, and K. G. Schwab, Nature 443, 
193 (2006). 

^ A. A. Clerk, F. Marquardt, and K. Jacobs, New J. Phys. 

10, 095010 (2008). 
® A. D. Armour and M. P. Blencowe, New J. Phys. 10, 

095004 (2008). 

M. D. LaHaye, J. Suh, P. M. Echternach, K. G. Schwab, 
and M. L. Roukes, Nature 459, 960 (2009). 

* M. Blencowe, Phys. Rep. 395, 159 (2004). 

® D. A. Rodrigues, J. Imbers, and A. D. Armour, Phys. Rev. 

Lett. 98, 067204 (2007). 
^° D. A. Rodrigues and A. D. Armour, New. J. Phys. 7, 251 
(2005). 

" A. D. Armour, M. P. Blencowe, and Y. Zhang, Phys. Rev. 
B 69, 125313 (2004). 

^2 J. Koch and F. von Oppen, Phys. Rev. Lett. 94, 206804 
(2005); J. Koch, M. E. Raikh, and F. von Oppen, 
Phys. Rev. Lett. 95, 056801 (2005); J. Koch, F. von Op- 
pen, and A. V. Andreev, Phys. Rev. B 74, 205438 (2006). 
D. Fedorets, L. Y. Gorelik, R. I. Shekhter, and M. Jonson, 
Phys. Rev. Lett. 92, 166801 (2004). 

" A. D. Armour and A. MacKinnon, Phys. Rev. B 66, 035333 
(2002). 

T. Novotny, A. Donarini, and A. -P. Jauho, 
Phys. Rev. Lett. 90, 256801 (2003). 

T. Novotny, A. Donarini, C. Flindt, and A.-P. Jauho, 

Phys. Rev. Lett. 92, 248302 (2004). 

T. Brandes, Phys. Rep. 408, 315 (2005). 



A. Braggio, J. Konig, and R. Fazio, Phys. Rev. Lett. 
96, 026805 (2006); C. Flindt, T. Novotny, A. Braggio, 
M. Sassetti, and A.-P. Jauho, Phys. Rev. Lett. 100, 150601 
(2008). 

A. Mitra, I. Aleiner, and A. J. Millis, Phys. Rev. B 69, 
245302 (2004). 

^° L. Glazman and R. Shehkter, Sov. Phys. JETP 67, 163 
(1988). 

N. S. Wingreen, K. W. Jacobsen, and J. W. Wilkins, 

Phys. Rev. B 40, 11834 (1989). 
^2 K. Flensberg, Phys. Rev. B 68, 205323 (2003). 
^3 A. A. Glerk, Phys. Rev. B 70, 245306 (2004). 

D. Mozyrsky, M. B. Hastings, and I. Martin, Phys. Rev. B 

73, 035104 (2006). 

R. P. Feynman and F. L. Vernon, Annals of Physics 24, 
118 (1963). 

^"^ D. Mozyrsky, I. Martin, and M. B. Hastings, Phys. Rev. 

Lett. 92, 018303 (2004). 
2^ S. D. Bennett and A. A. Glerk, Phys. Rev. B 74, 201301 

(2006). 

L. S. Schulman, Techniques and Applications of Path In- 
tegration (Dover Publications, Inc., 2005). 
U. Weiss, Quantum Dissipative Systems, vol. 13 (World 
Scientific Publishing, 2008). 
^° M. A. Armen and H. Mabuchi, Phys. Rev. A 73, 063801 
(2006). 

D. A. Rodrigues, J. Imbers, T. Harvey, and A. D. Armour, 
New. J. Phys. 9, 84 (2007). 

S. Strogatz, Nonlinear Dynamics and Chaos: With Ap- 
plications to Physics, Biology, Chemistry and Engineering 
(Westview Press, 2000). 

N. G. Van Kampen, Stochastic processes in physics and 
chemistry (Elsevier, 2008), 3rd ed. 

H. J. Haug and A.-P. Jauho, Quantum Kinetics m 
Transport and Optics of Semiconductors (Springer- Verlag, 
2008), 2nd ed. 



